%% fftGeneral
% 
% function to take fft of input signal and output fourier harmonics or
% amplitude, and phase.
%
%% Synatax
% [outputs] = fftGeneral(PressureData, Time, Plot_Type)
%
%
%
%% Description
% This function will take in data from pressure and flow extraction/fitting
% and run the fouier decompostion. The output of this function will be the
% fourier harmonics and amplitude and phase. Set plot_type to string of
% 'Pressure' for the pressure fft output and 'Flow' for the flow fft
%  output.

%% Arguments
%* PressureData                         - Double: A nx1 double of the
%                                         pressure data
%
%* Time                                 - Double: A nx1 double of time data
%
%
%* plot_type                            - String: A string of pressure or
%                                         flow to indicate type of fft plots to genrate
%
%% Returns
%
%* Outputs                              - Double: A nx5 double of
%                                         freqeuncy, freq(m) =  frequency components
%                                         a0 = DC offset
%                                         an = amplitude of the consines (real part, rec corrdinates)
%                                         bn = amplitude of the sines (imaginary part, rec corrdinates)
%                                         amp = amp (polar)
%                                         phase_m = phase (polar)
%
%% Function Side Effects
%* N/A
%
%% Exceptions
%* An error is thrown if data type of input varaibles is not compatible
%  with what the function is expecting
%
%% See Also
% N/A
 


function [outputs, fig1, t, ydata, z1, y, Yfft, delta_t, L, Tp] = fftGeneralWithMoreOutputs(ydata, t, plot_type)

if isa(ydata, 'double') ~=1
    type = class(ydata);
    error(sprintf(('Pressure Data is of type %s please change to type double'), type))    
end

if isa(t, 'double') ~=1
    type = class(t);
    error(sprintf(('Time Values are of type %s please change to type double'), type))
end
    
if isa(plot_type, 'char') ~=1
    type = class(plot_type);
    error(sprintf(('Plot Type selection is of of type %s please change to type char'), type))
end

%Set Figure Display Setttings
set(0, 'DefaultFigureVisible', 'off');




%% STEP 1
%------------Fourier decomposition of the wave analysis------------%   
%----------------this is a standard method-------------------------%



Y=ydata; % Column of Data Values
L=length(Y);
N=floor((L)/2);
delta_t=t(2)-t(1); % time step
Tp=t(end)+delta_t; % Time period is t(n+1)
Fs=1/Tp; % Sampling Frequency
Yfft = fft(Y);
f = [1:1:L-1]'./Tp;
NN=N-1;
Y1=Yfft(1:(L-1));
Y=Y1/L; %scale by n
Tp1=Tp;
z1=0:Tp1/L:(L-1)*Tp1/L;
i=sqrt(-1);
om=2*pi/Tp1;


%% STEP 2
%Fourier coefficients (cartesian complex plane an and bn) 
%standard method, shouldn't really need to be changed or adapted
a0=Y(1);
an=2.*real(Y(2:N));bn=-2.*imag(Y(2:N));

%Fourier coefficients in their polar coordinates 
amp=sqrt(an.^2 + bn.^2);
phase_m = atan(bn./an);
phase_m = unwrap(phase_m);

%phase_m = unwrap(phase_m); %This rescales the phase shift so that there are no unreasonable spikes
%This way of calculating the phase means that sometimes it can flick from side to side (I.e go from 2pi/3 to ?pi/3 between harmonics, so need to make sure it is smooth

fig1 = figure
plot((1:N-1),phase_m)



%This will then recomponse the wave, but also calculate the frequency of each harmonic? it is just m*omega /(2*pi)
%where m is the harmonic 
y=a0;
for m=1:NN
    Z_m=real((an(m)+i.*bn(m)).*exp(-i.*m.*om.*z1))';
    y=y+Z_m;
    freq(m)=m.*om/(2*pi);
end

%% STEP 3
%Plot the recomponsition of the wave to ensure that the decomposition is correct
%good check but not always necessary

if strcmp(plot_type,'Pressure')
    plot(t,ydata,'-b','LineWidth',2)
    hold on
    plot(z1,y,'-*r','LineWidth',2)
    xlabel('Time');ylabel('Y data Magnitude (dynes/cm^2')
    hold on
    inv_Y=ifftn(Yfft);
    ninv=L-1;Tp=Tp -delta_t;
    plot(0:Tp/ninv:(ninv-1)*Tp/ninv,inv_Y(1:floor(L-1)),'*k')
    title('Input Pressure Wave: Fourier Expansion')
else
    plot(t,ydata,'-b','LineWidth',2)
    hold on
    plot(z1,y,'-*r','LineWidth',2)
    xlabel('Time');ylabel('Y data Magnitude (cm/s)')
    hold on
    inv_Y=ifftn(Yfft);
    ninv=L-1;Tp=Tp -delta_t;
    plot(0:Tp/ninv:(ninv-1)*Tp/ninv,inv_Y(1:floor(L-1)),'*k')
    title('Input Flow Wave: Fourier Expansion')
end

legend('Data Points','Fourier Expansion','Inverse Fourier')


freq_new = [0, freq];
amp_new = [a0; amp];
phase_m_new = [0; phase_m];
an_new = [0; an];
bn_new = [0; bn];

set(0, 'DefaultFigureVisible', 'on');
outputs = [freq_new',amp_new,phase_m_new,an_new,bn_new];